%Plots for draft of TSI reconstruction manuscript. Constructed after
%revisions to the Gibbs sampler in accordance with the ACRIM Gap manuscript
%
% Ted Amdur
% 10/25/22 

clearvars

tsiComparison = 1; %Set to 1 to plot Figure 1 of manuscript
cycleMin = 0; %Set to 1 to plot Figure 2 of manuscript
tsiSunspots=0; %Set to 1 to plot Figure 3 of manuscript
tsiSunspotsDetail=0; %Set to 1 to plot detail of Figure 3 of manuscript
cycleShift=0; %Set to 1 to plot change in TSI between cycle minima
tsiComparisonPanelC = 0; %Set to 1 to plot only panel c of manuscript

fSize = 8;

%from runchain_22_04_25.m

load('ar2_24_12_10_long.mat')
%load('ar2_23_05_10_nrltsicomp.mat')
%load('ar2_23_05_10_pmodcorrections.mats')
load(outDat.obsmatrix); %From makeobsmatrix.m
excludeFliers=outDat.opts.excludeFliers;
dateStruct = getdates;
startDate = datejd(dateStruct.all(1));
endDate = datejd(dateStruct.all(2));
dates = [juliandate(datetime(1980,2,1)) juliandate(datetime(2023,12,31))]; %Use full interval
%Load data, with colLabels corresponding to observer source for each column
valMAll = valM;dateMAll=dateM;
if excludeFliers %excludeFliers
    load excludeMask_24_12_10.mat
    valM(excludeMask) = NaN;
    oM(excludeMask) = false;
    dateM=dateMAll;
end

%Reorder to match labels
%Order for proxies is sunspots first then MgII
%Order for sats is chronological from first observation: HF, ACRIM1, ERBE, ACRIM2,
%VIRGO/SOHO, ACRIM3, TIM/SORCE, PREMOS/PICARD, TCTE, TSIS-1
%Order as proxies followed by satellites in chronological order
lI=[8;4;6;1;5;2;12;3;9;7;10;11];
%sO=[7;3;6;4;12;5;9;8;10;11];
%Reorder arrays used in this plotting script
A=A(lI,:,:);
%Create alternative column labels for publication
colLabels=[...
    "SILSO";
    "Bremen Mg-II";
    "Nimbus-7/HF";
    "SMM/ACRIM1";
    "ERBS/ERBE";
    "UARS/ACRIM2";
    "SOHO/VIRGO"
    "ACRIMSAT/ACRIM3";
    "SORCE/TIM";
    "Picard/PREMOS";
    "TCTE/TIM";
    "TSIS-1/TIM";
    ];
offsets=offsets(lI);
oM=oM(:,lI);
sigY=sigY(lI,:);
tau=tau(:,lI);
valM=valM(:,lI); valMAll=valMAll(:,lI);
scaling=outDat.scaling(lI);
pindex=~outDat.satindex(lI);
oindex=outDat.oindex(lI);


%Colormap for plots
c1 = [51; 136; 68; 17; 153; 221; 204; 136; 170];
c2 = [34; 204; 170; 119; 153; 204; 102; 34; 68];
c3 = [136; 238; 153; 51; 51; 119; 119; 85; 153];
c = [c1 c2 c3]; c = c./255;
c(10,:) = c(1,:);

%Backups
cc = colororder;

if tsiComparison 
showTrend = 1;
smoothWindow = 6; %set smoothing (months)
%Create x-axis points for cycle demarcation
c21=datejd([dateStruct.cycles(1,:), fliplr(dateStruct.cycles(1,:))]);
c23=datejd([dateStruct.cycles(3,:), fliplr(dateStruct.cycles(3,:))]);
c25=datejd([dateStruct.cycles(5,:), fliplr(dateStruct.cycles(5,:))]);
yB = 450;
width=17.8;
height=13;
f = figure2('Position',[1 1 width.*50 height.*50]);
set(gcf, 'Color', 'w');
subplot('position',[.11 .78 .795 .18]) %Plot of proxy observations
yyaxis left
fill(c21,[0 0 yB yB],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c23,[0 0 yB yB],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c25,[0 0 yB yB],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
h(1)=plot(dateM,valM(:,1)+offsets(1),'.','MarkerSize',6);
ylabel('Sunspot number')
yyaxis right
h(2)=plot(dateM,valM(:,2)+offsets(2),'.','MarkerSize',6);
ylabel('Mg-II index')
xlim([datetime(1978,1,1) datetime(2024,1,1)])
yyaxis left
ylim([0 yB])
%Add labels of cycles and cycle minima
text(datejd(dateStruct.cycles(1,1))+years(3.5),yB-50,'Cycle 21','FontSize',fSize)
text(datejd(dateStruct.cycles(2,1))+years(2.25),yB-50,'Cycle 22','FontSize',fSize)
text(datejd(dateStruct.cycles(3,1))+years(3.75),yB-50,'Cycle 23','FontSize',fSize)
text(datejd(dateStruct.cycles(4,1))+years(3.25),yB-50,'Cycle 24','FontSize',fSize)
%Get low sunspot months and plot as lines:
[~,~,~,~,lowI]=getlowspotindex(dateM,valM,offsets,colLabels,5);
for i = find(lowI)
    % Extract current month
    currDate = dateM(i);
    
    % Calculate start and end of month
    startDate = dateshift(currDate, 'start', 'month');
    endDate = dateshift(currDate, 'end', 'month');
    
    % Overlay line segments on existing plot
    
    hold on;
    plot([startDate, endDate], [yB-50, yB-50], 'k-', 'LineWidth', 10);
    hold off;
end

set(gca,'FontSize',fSize)
text(datetime(1977,12,1),yB+25,'A','FontSize',fSize)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot of satellites
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot('position',[.11 .515 .795 .225]) %Plot of satellite observations
fill(c21,[1360 1360 1374 1374],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c23,[1360 1360 1374 1374],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c25,[1360 1360 1374 1374],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
ind = 1;
for ii = 3:numel(lI) %Iterate over satellite observations
    if showTrend
        hold on
        [trends,offsets2]= returntrend(A,tau,ii);
        tM = mean(trends,2) + mean(offsets2)+offsets(ii);
        t995 = quantile(trends,.975,2)+quantile(offsets2,.975)+offsets(ii);
        t005 = quantile(trends,.025,2)+quantile(offsets2,.025)+offsets(ii);
        x2 = [dateM(oM(:,ii))', flip(dateM(oM(:,ii)))'];
        fill(x2,[t995(oM(:,ii))',flip(t005(oM(:,ii)))'], ...
            [1 .85 .85],'FaceAlpha',0.5,'LineStyle','none');
        hold on
        plot(dateM(oM(:,ii)),tM(oM(:,ii)),'Color','r')
        hold on
    end
    plot(dateMAll,valMAll(:,ii)+offsets(ii),'o','MarkerSize',2.2,...
        'Color',c(ind,:));
    hold on
    hh(ind) = plot(dateM,valM(:,ii)+offsets(ii),'.','MarkerSize',6,...
        'Color',c(ind,:));
    ind = ind + 1;
end
xlim([datetime(1978,1,1) datetime(2024,1,1)])
ylim([1360 1374])
legend(hh,colLabels(3:end),'NumColumns',2)
legend boxoff
ylabel('TSI (W/m^{2})')
set(gca,'FontSize',fSize)
text(datetime(1977,11,1),1374.7,'B','FontSize',fSize)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot of reconstructions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot('position',[.11 .06 .795 .407])
%Set bounds of figure so other lines can riff on it
xB = [datetime(1978,1,1) datetime(2024,1,1)];
yB = [-0.35 2.25];
fill(c21,[yB(1)-.1 yB(1)-.1 yB(2)+.1 yB(2)+.1],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c23,[yB(1)-.1 yB(1)-.1 yB(2)+.1 yB(2)+.1],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c25,[yB(1)-.1 yB(1)-.1 yB(2)+.1 yB(2)+.1],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
cColor = get(gca,'colororder');
%Get CI of our estimate,plot
tsix = prctile(xAll',[2.5 10 50 90 97.5])';
for iS = 1:size(tsix,2)
    tsix(:,iS) = smoothPH(tsix(:,iS),smoothWindow);
end
[~,tsi0] = meanmin(dateM,tsix(:,3),23);
tsix = tsix - tsi0;
x2 = [dateM', fliplr(dateM')];
fill(x2,[tsix(:,1)',fliplr(tsix(:,end)')], [.85 .85 .85],'FaceAlpha',...
    0.5,'LineStyle','none');
hold on
fill(x2,[tsix(:,2)',fliplr(tsix(:,4)')], [.75 .75 .75],'FaceAlpha',0.5,...
    'LineStyle','none');

%     %Get satellite only reconstruction
load('ar2_24_12_12_satonly.mat')
%load('ar2_23_05_10_satonlydrift.mat')
%load('ar2_22_11_04_satnodrift.mat')
xms = mean(xAll,2);
xms = smoothPH(xms,smoothWindow);
xms = meanmin(dateM,xms,23);

% %Get MLR model reconstruction
if ~exist('xMLR','var')
    [xMLR,~] = makexmlr(dateM,pindex,oindex,valM,xAll,offsets);
end
% load('ar2_23_05_10_nrltsicomp.mat')
xmmlr=mean(xMLR,2);
xmmlr=smoothPH(xmmlr,smoothWindow);
% [~,xmmlrO] = meaninterval(dateMAll,xmmlr,1990,2010);
% xmmlr = xmmlr-xmmlrO;
[xmmlr] = meanmin(dateMAll,xmmlr,23);

%Plot other TSI reconstructions
ind = 1;
load oTSI_24_11_03.mat
load('lowspots_cutoff5_24_12_9.mat') 
%Plot PMOD
hold on
tsiA = smoothPH(oTSI(7).mthtsi,smoothWindow);
% [~,tsiAO] = meaninterval(oTSI(7).mthdatetime,tsiA,1990,2010);
tsiA0 = mean(tsiA(ismember(oTSI(7).mthdatetime,dateM(low23))));
h(ind) = plot(oTSI(7).mthdatetime,tsiA-tsiA0,'LineWidth',1.0,...
    'Color',c(7,:));
legendtxt(ind) = string(oTSI(7).product);
ind = ind + 1;

%Plot SOLID
hold on
tsiA = smoothPH(oTSI(9).mthtsi,smoothWindow);
% [~,tsiAO] = meaninterval(oTSI(9).mthdatetime,tsiA,1990,2010);
tsiA0 = mean(tsiA(ismember(oTSI(9).mthdatetime,dateM(low23))));
h(ind) = plot(oTSI(9).mthdatetime,tsiA-tsiA0,'LineWidth',1.0,...
    'Color',c(8,:));
legendtxt(ind) = "Comm.-Consensus";
ind = ind + 1;

%Plot NRLTSI2
hold on
tsiA = smoothPH(oTSI(4).mthtsi,smoothWindow);
% [~,tsiAO] = meaninterval(oTSI(4).mthdatetime,tsiA,1990,2010);
tsiA0 = mean(tsiA(ismember(oTSI(4).mthdatetime+hours(12),dateM(low23))));
h(ind) = plot(oTSI(4).mthdatetime,tsiA-tsiA0,'LineWidth',1.0,...
    'Color',c(5,:));
legendtxt(ind) = string(oTSI(4).product);
ind = ind + 1;


%Plot SATIRE-S during whole interval
hold on
tsiA = smoothPH(oTSI(5).mthtsi,smoothWindow);
% [~,tsiAO] = meaninterval(oTSI(5).mthdatetime,tsiA,1990,2010);
tsiA0 = mean(tsiA(ismember(oTSI(5).mthdatetime,dateM(low23))));
dateMDIHMI=oTSI(5).mthdatetime(oTSI(5).mthdatetime>datetime(1999,02,02));
%Plot SATIRE-S during period of high-quality MDI/HMI data
h(ind) = plot(dateMDIHMI,tsiA(oTSI(5).mthdatetime>datetime(1999,02,02))-tsiA0,'LineWidth',1.0,...
    'Color',c(4,:));
legendtxt(ind)=string(oTSI(5).product);
%Plot SATIRE-S during regression-dominated period
ind=ind+1;
h(ind)=plot(oTSI(5).mthdatetime,tsiA-tsiA0,'--','LineWidth',1.0,...
    'Color',c(4,:));
legendtxt(ind)="SATIRE-S regression";
hold on 

%Plot CMIP7 dataset
hold on
tsiA = smoothPH(oTSI(11).mthtsi,smoothWindow);
tsiA0 = mean(tsiA(ismember(oTSI(11).mthdatetime,dateM(low23))));
h(ind) = plot(oTSI(11).mthdatetime,tsiA-tsiA0,'LineWidth',1.0,...
    'Color',cc(2,:));
legendtxt(ind) = "CMIP7";

ind = ind + 1; %Index our contribution
hold on
indMain = ind;
h(ind) = plot(dateM,tsix(:,3),'LineWidth',1.75,'Color','k');
legendtxt(ind) = "BTSI";
ind = ind + 1;
hold on
h(ind)=plot(dateMAll,xms,'LineWidth',1.0,'Color',c(1,:));
legendtxt(ind)="BTSI sat. only";
ind=ind+1;
hold on
h(ind)=plot(dateMAll,xmmlr,'LineWidth',1.0,'Color',c(2,:));
hold on
h(indMain) = plot(dateM,tsix(:,3),'LineWidth',1.75,'Color','k');
line([dateMAll(1), dateMAll(end)], [0, 0], 'LineStyle', ':', 'Color', 'k')
legendtxt(ind)="BTSI proxy model";
legend(h,legendtxt,'NumColumns',3)
legend boxoff
set(gca,'FontSize',fSize)
xlabel('Year')
ylabel('TSI anomaly from Cyc. 23 minimum (W/m^{2})')
xlim(xB)
ylim(yB)
text(xB(1)-day(60),yB(2)+0.06,'C','FontSize',fSize)
[~,~,~,pthDate]=datechars;
    savePth=['plots/tsicompare_' pthDate '.eps'];
%export_fig savePth
x0=1;
y0=1;
set(gcf,'units','centimeters','position',[x0,y0,width,height])

exportgraphics(f,savePth,'ContentType','vector',...
               'BackgroundColor','none')
% exportgraphics(f,savePth,'ContentType','vector',...
%                'BackgroundColor','none')
end
if cycleMin
    %First, calculate the percentage of time that each simulation spends
    %inside of the 95% confidence interval

    xAll=xAll'+offsets(9);xms=mean(xAll,1);
    width=7.4;
    height=8.4;
    f = figure2('Position',[100 -242 width.*100 height.*100]);
    set(gcf, 'Color', 'w');
    subplot('position',[.16 .42 .8 .54])
    trendInd=dateM>=datejd(dates(1))&dateM<datejd(dates(2));
    X = [ones(size(dateM(trendInd),1),1) juliandate(dateM(trendInd))];
    warning('off','MATLAB:rankDeficientMatrix') %Some realizations for quant_reg are rank deficient
    for ii=1:size(xAll,1)
        b5(:,ii)=quant_reg(X,xAll(ii,trendInd)',0.05);
        b50(:,ii)=quant_reg(X,xAll(ii,trendInd)',0.50);
        yhat5(:,ii) = X*b5(:,ii);
        yhat50(:,ii) = X*b50(:,ii);
    end
    hold on
    yplot5=prctile(yhat5',[2.5 10 50 90 97.5])';
    yplot50=prctile(yhat50',[2.5 10 50 90 97.5])';
    x2 = [dateM(trendInd)', fliplr(dateM(trendInd)')];
    fill(x2,[yplot5(:,1)',fliplr(yplot5(:,5)')],c(2,:).*1.05,'FaceAlpha',...
        0.5,'LineStyle','none');
    fill(x2,[yplot5(:,2)',fliplr(yplot5(:,4)')],c(2,:).*.9,'FaceAlpha',...
        0.5,'LineStyle','none');
    hold on
    plot(dateM(trendInd),yplot5(:,3),'Color',c(2,:).*.9,'LineWidth',2)
    hold on
    fill(x2,[yplot50(:,1)',fliplr(yplot50(:,5)')],c(8,:).*1.05,'FaceAlpha',...
        0.5,'LineStyle','none');
    fill(x2,[yplot50(:,2)',fliplr(yplot50(:,4)')],c(8,:).*.9,'FaceAlpha',...
        0.5,'LineStyle','none');
    hold on
    plot(dateM(trendInd),yplot50(:,3),'Color',c(8,:),'LineWidth',2)
    hold on
    plot(dateM,xAll(1:100:end,:)','.','MarkerSize', 2, ...
        'Color',[0.65 0.65 0.65])
    hold on
    plot(dateM,xms,'Color','k','LineWidth',1.25)
    xlabel('Year')
    ylabel('TSI (W/m^{2})')
    set(gca,'FontSize',fSize)
    box off
    yticks(1360:1:1363)
    hA = get(gca,'YAxis');
    set(gca,'YMinorTick','on')
    hA.MinorTick = 'on';
    hA.MinorTickValues = 1360:0.25:1363;
    ylim([1360 1363])
    text(datetime(1979,1,1),1363.065,'A','FontSize',fSize)
    
    %Plot histograms
    subplot('position',[.16 .11 .8 .22])
    histogram(b5(2,:).*(365.25*10),linspace(-.40,0.02,85), ...
        'FaceColor',c(2,:), 'EdgeColor', 'w')
    hold on
    histogram(b50(2,:).*(365.25*10),linspace(-.40,0.02,85), ...
        'FaceColor',c(8,:), 'EdgeColor', 'w')
    ylabel('Realizations')
    xlabel('W/m^{2}/decade')
    xlim([-0.25 0])
    set(gca,'FontSize',fSize)
    text(-0.249,size(xAll,1)./4.2,'B','FontSize',fSize)
    [~,~,~,pthDate]=datechars;
    savePth=['plots/mincycle_' pthDate '.eps'];
    %export_fig savePth
    x0=1;
    y0=1;
    width=8.7;
    height=9.5;
    set(gcf,'units','centimeters','position',[x0,y0,width,height])

    exportgraphics(f,savePth,'ContentType','vector',...
        'BackgroundColor','none')
    %export_fig savePth
end
if tsiSunspots
    xAll=xAll-nanmean(xAll(:));
    xm = nanmean(xAll,2);
    yM=valM(:,1)+offsets(1);
    nanI=isnan(yM);
    y=yM(~nanI);x=xm(~nanI);
    %Determine piecewise regression for tsi-sunspot relationship
    %Code derived from following post:
    %https://www.mathworks.com/matlabcentral/answers/344059-piece-wise-linear-fitting
    %Assume the curve is defined in terms of two linear segments, with break at b1
    plusfun = @(x) max(x,0);
    model = @(P,x) P(1) + P(2)*x + (P(3)-P(2))*plusfun(x-P(4));
    %Parameter set P includes: offset, slope 1, slope 2, break value
    lb=[-50,-500,-500,min(x)+1E-8];
    ub=[50,500,500,max(x)-1E-8];
    P0=(lb+ub)./2; %naive, unbiased initial search location
    [Pfit,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(model,P0,x,y,lb,ub);
    conf = nlparci(Pfit,residual,'jacobian',jacobian);
    [Ypred,delta] = nlpredci(model,x,Pfit,residual,'Jacobian',jacobian);
    modelpred=model(Pfit,x);
    lowpred=modelpred+delta;
    highpred=modelpred-delta;
    [linfit,linint]=regress(y,[ones(length(x),1) x]);
    [~,sI]=sort(x); %Plot line in order on the x-axis
    
    f = figure2('Position',[10 10 1000 1000]);
    set(gcf, 'Color', 'w');
    x2=[x(sI); flipud(x(sI))];
    
    %Plot CI from BTSI for sunspots
    spotPred=zeros(size(x,1),size(xAll,2))';
    for iR=1:size(xAll,2)
        spotPred(iR,:)=squeeze(A(1,[2 3],iR))*[x';tau(~nanI,1)'].*scaling(1)+offsets(1);
    end
    linPreds=prctile(spotPred,[2.5 50 97.5])';
    linlow=linPreds(:,1);
    linhigh=linPreds(:,3);
    
    fill(x2,[linlow(sI); flipud(linhigh(sI))], [.85 .85 .85],'FaceAlpha',...
        0.5,'LineStyle','none');
    hold on
    fill(x2,[lowpred(sI); flipud(highpred(sI))],[1 .85 .85],'FaceAlpha',...
        0.5,'LineStyle','none');
    hold on
    
    
    
    cColor = colororder;
    for ii = 1:size(dateStruct.cycles,1)
        cycleI = dateM > datejd(dateStruct.cycles(ii,1)) & dateM < datejd(dateStruct.cycles(ii,2));
        hold on
        h(ii) = plot(xm(cycleI),valM(cycleI,1)+offsets(1),'.','MarkerSize',8,'Color',c(2*ii-1,:));
    end
    hold on
    h(ii+1)=plot(x(sI),linPreds(sI,2),'Color','k','LineWidth',1.5);
    hold on
    h(ii+2)=plot(x(sI),modelpred(sI),'Color','r','LineWidth',1.5);
    hold on
    rectangle('Position', [-.754, 0, 0.5, 50], ...
                'FaceColor', 'none', ...
                'EdgeColor', 'k');
    legend(h,'Cycle 21','Cycle 22','Cycle 23','Cycle 24','Cycle 25','BTSI prediction','Piecewise prediction','Location','SouthEast')
    legend boxoff
    hold on
    xlim([-1 1.35])
    ylim([0 300])
    xlabel('BTSI anomaly (W/m^{2})')
    ylabel('SILSO sunspot number')
    set(gca,'FontSize',fSize)
    [~,~,~,pthDate]=datechars;
    savePth=['plots/tsispots_' pthDate '.eps'];

    x0=1;
    y0=1;
    width=17.8;
    height=17.8;
    set(gcf,'units','centimeters','position',[x0,y0,width,height])

    exportgraphics(f,savePth,'ContentType','vector',...
        'BackgroundColor','none')
    
end
if tsiSunspotsDetail
    xAll=xAll-nanmean(xAll(:));
    xm = nanmean(xAll,2);
    yM=valM(:,1)+offsets(1);
    nanI=isnan(yM);
    y=yM(~nanI);x=xm(~nanI);
    %Determine piecewise regression for tsi-sunspot relationship
    %Code derived from following post:
    %https://www.mathworks.com/matlabcentral/answers/344059-piece-wise-linear-fitting
    %Assume the curve is defined in terms of two linear segments, with break at b1
    plusfun = @(x) max(x,0);
    model = @(P,x) P(1) + P(2)*x + (P(3)-P(2))*plusfun(x-P(4));
    %Parameter set P includes: offset, slope 1, slope 2, break value
    lb=[-50,-500,-500,min(x)+1E-8];
    ub=[50,500,500,max(x)-1E-8];
    P0=(lb+ub)./2; %naive, unbiased initial search location
    [Pfit,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(model,P0,x,y,lb,ub);
    conf = nlparci(Pfit,residual,'jacobian',jacobian);
    [Ypred,delta] = nlpredci(model,x,Pfit,residual,'Jacobian',jacobian);
    modelpred=model(Pfit,x);
    lowpred=modelpred+delta;
    highpred=modelpred-delta;
    [linfit,linint]=regress(y,[ones(length(x),1) x]);
    [~,sI]=sort(x); %Plot line in order on the x-axis
    
    f = figure2('Position',[10 10 1000 1000]);
    set(gcf, 'Color', 'w');
    x2=[x(sI); flipud(x(sI))];
    
    %Plot CI from BTSI for sunspots
    spotPred=zeros(size(x,1),size(xAll,2))';
    for iR=1:size(xAll,2)
        spotPred(iR,:)=squeeze(A(1,[2 3],iR))*[x';tau(~nanI,1)'].*scaling(1)+offsets(1);
    end
    linPreds=prctile(spotPred,[2.5 50 97.5])';
    linlow=linPreds(:,1);
    linhigh=linPreds(:,3);
    
    fill(x2,[linlow(sI); flipud(linhigh(sI))], [.85 .85 .85],'FaceAlpha',...
        0.5,'LineStyle','none');
    hold on
    fill(x2,[lowpred(sI); flipud(highpred(sI))],[1 .85 .85],'FaceAlpha',...
        0.5,'LineStyle','none');
    hold on
    
    
    
    cColor = colororder;
    for ii = 1:size(dateStruct.cycles,1)
        cycleI = dateM > datejd(dateStruct.cycles(ii,1)) & dateM < datejd(dateStruct.cycles(ii,2));
        hold on
        h(ii) = plot(xm(cycleI),valM(cycleI,1)+offsets(1),'.','MarkerSize',8,'Color',c(2*ii-1,:));
    end
    hold on
    h(ii+1)=plot(x(sI),linPreds(sI,2),'Color','k','LineWidth',1.5);
    hold on
    h(ii+2)=plot(x(sI),modelpred(sI),'Color','r','LineWidth',1.5);
    set(gca,'xtick',[])
    set(gca,'ytick',[])
    xlim([-.754 -.754+.5])
    ylim([0 50])
    set(gca,'FontSize',fSize)
    [~,~,~,pthDate]=datechars;
    savePth=['plots/tsispotsdetail_' pthDate '.eps'];

    x0=1;
    y0=1;
    width=6.5;
    height=5.2;
    set(gcf,'units','centimeters','position',[x0,y0,width,height])

    exportgraphics(f,savePth,'ContentType','vector',...
        'BackgroundColor','none')
end
if cycleShift
%     %first, select for datetimes when SILSO spot count is less than 5
%     sptI=find(strcmp(colLabels,"SILSO")); %Get column of SILSO observations
%     lowI=(valM(:,sptI)+offsets(sptI))<5;
%     dateS=getdates;
%     low24=and(lowI,dateM.Year<2025 & dateM.Year > 2015);
%     low23=and(lowI,dateM.Year<2015 & dateM.Year > 2003);
%     low22=and(lowI,dateM.Year<2000 & dateM.Year > 1990);
%     low21=and(lowI,dateM.Year<1990 & dateM.Year > 1980);
    
    [low21,low22,low23,low24]=getlowspotindex(dateM,valM,offsets,colLabels,5);
    
    %Control for slight differences in low sunspot regimes by subtracting
    %expected TSI difference predicted by sunspots
    xCorr=xAll-((repmat(valM(:,1),[1 size(A,3)])-min(valM(:,1)))./...
        outDat.scaling(8))./squeeze(A(1,2,:))';
    diff2423=mean(xCorr(low24,:),1)-mean(xCorr(low23,:),1);
    diff2422=mean(xCorr(low24,:),1)-mean(xCorr(low22,:),1);
    diff2421=mean(xCorr(low24,:),1)-mean(xCorr(low21,:),1);
    f = figure2('Position',[110 -200 800 500]);
    set(gcf, 'Color', 'w');
    %First do end of cycle 24 minus 23
    pEdges = linspace(-0.35,0.05,1000);
    [y1,x1] = histcounts(diff2423,pEdges);
    [x1,y1,~] = histtoplot(y1,x1,50);
    hold on
    plot(x1,y1,'Color',c(1,:),'LineWidth',2.5)
    
    %Then do end of cycle 24 minus 22
    pEdges = linspace(-0.35,0.05,1000);
    [y2,x2] = histcounts(diff2422,pEdges);
    [x2,y2,~] = histtoplot(y2,x2,150);
    hold on
    plot(x2,y2,'Color',c(2,:),'LineWidth',2.5)
    
    %Then do end of cycle 24 minus 21
    pEdges = linspace(-0.35,0.05,1000);
    [y3,x3] = histcounts(diff2421,pEdges);
    [x3,y3,~] = histtoplot(y3,x3,150);
    hold on
    plot(x3,y3,'Color',c(3,:),'LineWidth',2.5)
    xlabel('\Delta TSI (W/m^{2})')
    ylabel('PDF')
    lgd = legend('Cycle 24 - Cycle 23','Cycle 24 - Cycle 22',...
        'Cycle 24 - Cycle 21','Location','NorthWest');
    legend boxoff
%     fontsize(lgd,8,'points');
    set(gca,'FontSize',fSize)
    [~,~,~,pthDate]=datechars;
    savePth=['plots/mindiff_' pthDate '.eps'];
    x0=1;
    y0=1;
    width=8.7;
    height=6;
    set(gcf,'units','centimeters','position',[x0,y0,width,height])

    exportgraphics(f,savePth,'ContentType','vector',...
        'BackgroundColor','none')
    
    
    
%     histogram(diff2423,'BinWidth',0.01)
%     hold on
%     histogram(diff2422,'BinWidth',0.01)
%     hold on
%     histogram(diff2421,'BinWidth',0.01)
%     legend('23','22','21')
end
if tsiComparisonPanelC
showTrend = 1;
smoothWindow = 6; %set smoothing (months)
%Create x-axis points for cycle demarcation
c21=datejd([dateStruct.cycles(1,:), fliplr(dateStruct.cycles(1,:))]);
c23=datejd([dateStruct.cycles(3,:), fliplr(dateStruct.cycles(3,:))]);
c25=datejd([dateStruct.cycles(5,:), fliplr(dateStruct.cycles(5,:))]);

f = figure2('Position',[1 1 750 500]);
set(gcf, 'Color', 'w');
fill(c21,[-1.2 -1.2 1.8 1.8],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c23,[-1.2 -1.2 1.8 1.8],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
fill(c25,[-1.2 -1.2 1.8 1.8],[.96 .96 .863],'FaceAlpha',...
    0.4,'LineStyle','none');
hold on
cColor = get(gca,'colororder');
%Get CI of our estimate,plot
tsix = prctile(xAll',[2.5 10 50 90 97.5])';
for iS = 1:size(tsix,2)
    tsix(:,iS) = smoothPH(tsix(:,iS),smoothWindow);
end
[~,tsiO] = meaninterval(dateM,tsix(:,3),1990,2010);
tsix = tsix-tsiO;
x2 = [dateM', fliplr(dateM')];
fill(x2,[tsix(:,1)',fliplr(tsix(:,end)')], [.85 .85 .85],'FaceAlpha',...
    0.5,'LineStyle','none');
hold on
fill(x2,[tsix(:,2)',fliplr(tsix(:,4)')], [.75 .75 .75],'FaceAlpha',0.5,...
    'LineStyle','none');

%     %Get satellite only reconstruction
load('ar2_24_11_22_satonly.mat')
%load('ar2_23_05_10_satonlydrift.mat')
%load('ar2_22_11_04_satnodrift.mat')
xms = mean(xAll,2);
xms = smoothPH(xms,smoothWindow);
[~,xmsO] = meaninterval(dateMAll,xms,1990,2010);
xms = xms-xmsO;

% %Get MLR model reconstruction
if ~exist('xMLR','var')
    [xMLR,~] = makexmlr(dateM,pindex,oindex,valM,xAll,offsets);
end
% load('ar2_23_05_10_nrltsicomp.mat')
xmmlr=mean(xMLR,2);
xmmlr=smoothPH(xmmlr,smoothWindow);
[~,xmmlrO] = meaninterval(dateMAll,xmmlr,1990,2010);
xmmlr = xmmlr-xmmlrO;

%Plot other TSI reconstructions
ind = 1;
load oTSI_24_11_03.mat
%Plot PMOD
hold on
tsiA = smoothPH(oTSI(7).mthtsi,smoothWindow);
[~,tsiAO] = meaninterval(oTSI(7).mthdatetime,tsiA,1990,2010);
h(ind) = plot(oTSI(7).mthdatetime,tsiA-tsiAO,'LineWidth',2.5,...
    'Color',c(7,:));
legendtxt(ind) = string(oTSI(7).product);
ind = ind + 1;

%Plot SOLID
hold on
tsiA = smoothPH(oTSI(9).mthtsi,smoothWindow);
[~,tsiAO] = meaninterval(oTSI(9).mthdatetime,tsiA,1990,2010);
h(ind) = plot(oTSI(9).mthdatetime,tsiA-tsiAO,'LineWidth',2.5,...
    'Color',c(8,:));
legendtxt(ind) = "Comm.-Consensus";
ind = ind + 1;

%Plot NRLTSI2
hold on
tsiA = smoothPH(oTSI(4).mthtsi,smoothWindow);
[~,tsiAO] = meaninterval(oTSI(4).mthdatetime,tsiA,1990,2010);
h(ind) = plot(oTSI(4).mthdatetime,tsiA-tsiAO,'LineWidth',2.5,...
    'Color',c(5,:));
legendtxt(ind) = string(oTSI(4).product);
ind = ind + 1;


%Plot SATIRE-S during whole interval
hold on
tsiA = smoothPH(oTSI(5).mthtsi,smoothWindow);
[~,tsiAO] = meaninterval(oTSI(5).mthdatetime,tsiA,1990,2010);
dateMDIHMI=oTSI(5).mthdatetime(oTSI(5).mthdatetime>datetime(1999,02,02));
%Plot SATIRE-S during period of high-quality MDI/HMI data
h(ind) = plot(dateMDIHMI,tsiA(oTSI(5).mthdatetime>datetime(1999,02,02))-tsiAO,'LineWidth',2.5,...
    'Color',c(4,:));
legendtxt(ind)=string(oTSI(5).product);
%Plot SATIRE-S during regression-dominated period
ind=ind+1;
h(ind)=plot(oTSI(5).mthdatetime,tsiA-tsiAO,'--','LineWidth',2.5,...
    'Color',c(4,:));
legendtxt(ind)="SATIRE-S regression";
hold on 

ind = ind + 1; %Index our contribution
hold on
h(ind) = plot(dateM,tsix(:,3),'LineWidth',3.5,'Color','k');
legendtxt(ind) = "BTSI";
ind = ind + 1;
hold on
h(ind)=plot(dateMAll,xms,'LineWidth',2.5,'Color',c(1,:));
legendtxt(ind)="BTSI sat. only";
ind=ind+1;
hold on
h(ind)=plot(dateMAll,xmmlr,'LineWidth',2.5,'Color',c(2,:));

%Get low sunspot months and plot as lines:
[~,~,~,~,lowI]=getlowspotindex(dateM,valM,offsets,colLabels,5);
for i = find(lowI)
    % Extract current month
    currDate = dateM(i);
    
    % Calculate start and end of month
    startDate = dateshift(currDate, 'start', 'month');
    endDate = dateshift(currDate, 'end', 'month');
    
    % Overlay line segments on existing plot
    hold on;
    plot([startDate, endDate], [-0.85, -0.85], 'k-', 'LineWidth', 12);
    hold off;
end


legendtxt(ind)="BTSI proxy model";
legend(h,legendtxt,'NumColumns',2)
legend boxoff
set(gca,'FontSize',fSize)
xlabel('Year')
ylabel('TSI anomaly from 1990-2010 mean (W/m^{2})')
xlim([datetime(1978,1,1) datetime(2024,1,1)])
ylim([-0.9 1.65])
text(datejd(dateStruct.cycles(1,1))+years(3.5),yB(1)+0.1,'Cycle 21','FontSize',fSize)
text(datejd(dateStruct.cycles(2,1))+years(2.25),yB(1)+0.1,'Cycle 22','FontSize',fSize)
text(datejd(dateStruct.cycles(3,1))+years(3.75),yB(1)+0.1,'Cycle 23','FontSize',fSize)
text(datejd(dateStruct.cycles(4,1))+years(3.25),-0.8,'Cycle 24','FontSize',fSize)
[~,~,~,pthDate]=datechars;
    savePth=['plots/tsicompareC_' pthDate '.eps'];
%export_fig savePth
x0=1;
y0=1;
width=17.8;
height=11.867;
set(gcf,'units','centimeters','position',[x0,y0,width,height])

exportgraphics(f,savePth,'ContentType','vector',...
               'BackgroundColor','none')
end




